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A  Semiclassical  Reactive  Flux  Method  for  the  Calculation  of 
Condensed  Phase  Activated  Rate  Constants 

Diane  E.  Sagnella,  Jianshu  Cao,  and  Gregory  A.  Voth 

Department  of  Chemistry ,  University  of  Pennsylvania,  Philadelphia,  Pennsylvania  19104- 
6323,  USA 

ABSTRACT 

A  semiclassical  reactive  flux  algorithm  for  calculating  thermally  activated  rate  con¬ 
stants  is  presented  which  is  based  on  a  semiclassical  transition  state  theory  due  to  Chap¬ 
man,  Garrett,  and  Miller  [J.  Chem.  Phys.  63  (1975)  2710].  This  reactive  flux  technique, 
when  combined  with  the  semiclassical  TST,  enables  one  to  describe  dynamical  recross¬ 
ings  of  the  transition  state  on  the  same  footing  as  tunneling  effects.  Most  importantly, 
the  method  is  readily  applied  to  nonlinear  multidimensional  systems  over  a  wide  range  of 
temperatures.  It  will  be  shown  that  the  method  works  very  well  for  a  variety  of  existing 
models. 
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1.  Introduction 

It  is  highly  desirable  to  have  computational  methods  which  allow  for  the  efficient  and 
accurate  numerical  evaluation  of  activated  rate  constants  in  complex  systems.  It  is  per¬ 
haps  no  surprise,  therefore,  that  the  preeminent  theory  of  this  sort -classical  transition 
state  theory  (TST)— is  one  of  the  central  theoretical  tools  in  modern  chemistry  [1].  Yet, 
classical  transition  state  theory  is  lacking  on  at  least  two  accounts.  First,  the  fundamental 
assumption  of  classical  TST  is  based  on  the  notion  that  the  recrossings  of  trajectories 
through  the  transition  state  can  be  neglected.  Second,  the  classical  TST  formulation  ne¬ 
glects  quantum  mechanical  effects  such  as  tunneling  and  mode  quantization.  Fortunately, 
the  influence  of  classical  recrossings  can  be  accounted  for  by  analytic  approaches  [2]  or  with 
the  reactive  flux  correlation  function  methodology  [3],  while  quantum  effects  can  be  de¬ 
scribed  by  semiclassical  [4-9]  or  quantum  mechanical  [10,11]  TST-like  approaches  (QTST). 
One  of  the  latter  QTST  approaches,  path  integral  QTST  [10],  has  also  recently  been  com¬ 
bined  by  ansatz  with  a  reactive  flux-like  formulation  so  that  both  recrossing  and  quantum 
effects  in  activated  rate  processes  can  be  estimated  [12].  A  different  approach  within  the 
path  integral  QTST  is  based  on  an  analytic  theory  [13]  for  the  recrossing  corrections  (i.e., 
the  “transmission  coefficient”  problem),  but  the  full  numerical  implementation  of  the  path 
integral  QTST  expression  is  used  to  determine  the  quantum  corrections  to  the  classical 
TST  rate. 

In  the  present  paper,  a  different  semiclassical  TST  (SCTST)  formulation  proposed 
nearly  20  years  ago  by  Chapman,  Garrett,  and  Miller  (CGM)  [8]  is  revisited.  As  it  turns 
out,  it  is  quite  illuminating  to  re-examine  this  SCTST  formulation  in  the  light  of  several 
more  recent  developments  in  the  field  of  condensed  phase  activated  dynamics  [2,3].  To  be 
specific,  it  will  be  shown  that  the  CGM  semiclassical  theory  is  continuously  accurate  and 
stable  at  temperatures  both  above  and  below  the  so-called  “crossover”  temperature  (i.e., 
the  temperature  below  which  the  dynamics  becomes  dominated  by  tunneling).  Moreover, 
the  accuracy  of  the  CGM  formula  will  be  shown  to  be  improved  by  rotation  of  the  TST 
dividing  surface.  Finally,  and  most  importantly,  the  SCTST  formula  is  shown  to  form 
the  basis  for  a  semiclassical  reactive  flux  algorithm  [3]  which  allows  one  to  estimate  both 
the  quantum  and  dynamical  corrections  to  the  classical  TST  rate  constant  on  the  same 
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footing. 

The  present  paper  is  organized  as  follows:  In  Sec.  2,  the  basic  CGM  formulation 
of  a  SCTST  expression  is  reviewed  and  the  attributes  of  the  formula  are  discussed.  The 
semiclassical  reactive  flux  algorithm  based  on  the  SCTST  formula  will  also  be  described. 
Next,  in  Sec.  3  several  applications  of  the  method  will  be  presented  which  are  designed  to 
test  the  theory  against  known  results.  Concluding  remarks  are  then  given  in  Sec.  4. 

2.  Theory 

A.  Semiclassical  TST  of  Chapman,  Garrett,  and  Miller 

The  basis  for  the  SCTST  formulation  of  CGM  [8]  is  the  semiclassical  perspective 
on  Boltzmann  statistical  mechanics  developed  by  Miller  [14].  For  example,  the  diagonal 
matrix  elements  of  the  Boltzmann  density  matrix  (i.e.,  the  particle  density)  are  given  by 

p(q)  =  (q|e*p(-/?tf)|q>  (1) 

where  0  =  (fc^T)-1.  By  splitting  the  Boltzmann  operator  into  a  product  of  two  equal 
parts,  and  by  inserting  a  complete  set  of  states  between  those  parts,  the  particle  density 
can  be  re-written  as  [14] 


P( qo)  =  J dqi\(q0\e  ^W/2|qi>|2  •  (2) 

Upon  introduction  of  the  semiclassical  limit  for  the  imaginary  time  propagator  elements 
(qo|e_/3H,/2|qi),  the  particle  density  then  takes  the  form  [14] 

P(qo)  =  h~N  J  dqi\dqi/dpo\~l  expi^-^  driJ[p(r),  q(r)]  j  ,  (3) 

where  |dqi/<9po|  is  the  Jacobian  relating  the  coordinate  at  imaginary  time  t\  —  h0/2  to 
its  momentum  at  imaginary  time  to  =  0,  i?[p(r),q(r)]  is  the  value  of  the  Hamiltonian, 
and  N  is  the  dimensionality  of  the  system.  The  trajectories  which  are  used  to  calculate 
the  value  of  the  Hamiltonian  in  eq.  (3)  are  generated  by  the  differential  equations 
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q(r)  = 


P(r)  = 


dH 


(4) 


ap(r)  ’  ^  7  «q(r) 

subject  to  the  initial  conditions  q(0)  =  qo  and  p(0)  =  po-  From  eq.  (4),  it  is  seen  that 
the  trajectories  are  propagated  on  the  upside-down  potential,  -V(q),  and  thereby  have 
excursions  into  the  classically  forbidden  regions  of  the  potential. 

An  insightful  trick  due  to  Miller  [14]  is  to  perform  a  change  of  variables  from  the 
differential  final  position  elements  dqi  to  the  differental  initial  momentum  elements  dpo- 
The  Jacobian  factor  from  this  transformation  cancels  the  Jacobian  factor  in  eq.  (3)  and 
one  is  left  with  only  an  integration  over  initial  momenta.  This  step  in  itself  is  a  simplifying 
feature  since  the  calculation  of  the  semiclassical  canonical  density  no  longer  requires  the 
solution  of  a  classical  trajectory  problem  having  double  ended  boundary  conditions.  For 
example,  the  semiclassical  partition  function  in  Miller’s  theory  is  given  by  [14] 


Qsc  —  J  dqopsc(qo) 

=  h~N  J  J  dpodqoexp  j  -  jj-  drH [p(r),  q(r)]  j 


(5) 


Stratt  and  Miller  [15]  have  noted  that  the  underlying  semiclassical  phase  space  distri¬ 
bution  function  in  this  theory  is  given  by 


Wsc(Po,qo)  =  h  N exp  ^  - 


f  2  fh0/2 

{hi  ' 


drff[p(r),q(r)] 


} 


(6) 


Actually,  since  Liouville’s  theory  implies  the  conservation  of  phase  space  differential  ele¬ 
ment  [i.e.,  dqodpo  =  dq(r)dp(r)],  then  eq.  (6)  can  be  rewritten  for  an  arbitrary  time  shift 
tq  to  give  the  more  general  semiclassical  phase  space  distribution  function  [8] 


Wsc(po,qo)  =  h  N exp 


H 


fh/3/2—ro 

/  dTH[p(T),q(r)\ 

J  —TQ 


} 


(7) 


The  trajectories  used  to  calculate  i/[p(r),  q(r)]  in  this  distribution  function  begin  with 
values  (po,qo)  at  time  r  =  0,  then  are  propagated  on  the  upside-down  potential  [cf.  eq. 
(4)]  backward  in  time  to  r  =  —to  and  forward  in  time  to  t  =  ft/3  —  tq. 
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The  optimal  choice  of  to  in  eq.  (7)  will  depend  on  the  quantity  which  is  being  averaged. 
(The  semiclassical  partition  function  Qsc  in  eq-  (5),  however,  is  independent  of  the  choice 
of  to  [8].)  For  example,  if  one  calculates  the  semiclassical  average  ( q 2)  with  the  distribution 
in  eq.  (7)  for  a  simple  harmonic  oscillator  of  frequency  u/,  one  finds 


(<z2>  = 


ft,  sinh(h(3u  —  2u^ro)  +  sinh{2uTo) 


(8) 


2  mu  2sinh2(h(3uj/2) 

For  this  example,  only  the  choice  of  to  =  0  or  ftp/ 2  gives  the  correct  answer. 

In  general,  one  does  not  know  the  optimal  choice  of  to,  except  for  certain  analytically 
solvable  models.  Based  on  such  information,  however,  one  can  assume  a  certain  value  for 
To  when  using  the  semiclassical  distribution  function  for  non-trivial  systems.  It  is  in  this 
same  spirit  that  CGM  suggested  an  approximate  SCTST  expression  for  the  rate  constant 
[8].  That  expression  for  the  forward  rate  constant  of  a  general  system  is  given  by 


kscrsr  =  QRlh~N  /  ■  ■  ■  / dPdQdpdqHq'  -  g)~h(p) 


xexp 


(  2  Chp/ 


hp/2-ro 


drH\p{T),q(T),  P(t),  Q(t)] 


}■ 


(9) 


where  Qr  is  the  partition  function  of  the  system  in  the  reactant  state,  q  and  p  are,  respec¬ 
tively,  the  reaction  coordinate  and  its  conjugate  momentum,  q  =  q*  defines  the  dividing 
surface  along  the  reaction  coordinate,  h( ■  ■  ■)  is  the  unit  step  function,  and  Q  and  P  are  all 
other  coordinates  and  conjugate  momenta  of  the  system,  respectively.  The  subscript  “0” 
on  the  momenta  and  coordinates  in  the  integrals  have  been  dropped  here,  but  it  should  be 
remembered  that  these  variables  provide  the  initial  conditions  for  the  trajectories  which 
are  propagated  forward  and  backward  in  time  on  the  upside-down  potential  in  order  to 
calculate  H[p(T),q(T),  P(t),Q(t)]  in  eq.  (9). 

Based  on  their  analysis  of  the  simple  parabolic  barrier,  CGM  concluded  that  the  most 
reasonable  choice  of  to  in  eq.  (9)  is  To  =  h/3/ 4.  More  specifically,  this  choice  of  to  is  the 
only  one  that  gives  the  exact  rate  constant  for  the  parabolic  barrier  model.  Interestingly, 
this  result  for  the  parabolic  barrier  coincides  with  the  minimum  of  the  flux  in  eq.  (9)  with 
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respect  to  the  choice  of  tq.  We  therefore  suggest,  without  proof,  that  the  best  choice  of 
ro  in  the  SCTST  formula  [eq.  (9)]  is  the  one  that  minimimizes  the  imaginary  time  action 
in  the  exponential  of  eq.  (9).  For  example,  at  any  initial  point  in  the  reaction  coordinate 
phase  space  (p,  q ),  the  imaginary  time  action  for  a  one-dimensional  barrier  problem  is  given 
by 


/fi/3/2— r0  fhp/l-TQ 

H(t)cIt  =  E  +  2  /  dTVb[q(r)]  ,  (10) 

-To  J  —Tq 

where  E  is  the  conserved  energy  of  the  particle  moving  on  the  inverted  barrier  potential 
—Vb(q).  If  the  extremum  of  the  action  is  taken,  one  obtains  the  condition  for  ro: 

Hh(-ii)]  =  ntoWJ/2  -  id)]  ■  (ii) 

Thus,  the  condition  of  ro  =  h/3/A  can  be  derived  from  eq.  (11)  for  a  symmetric  barrier 
uncoupled  to  a  bath.  In  view  of  the  additional  complications  in  solving  eq.  (11)  for  the 
most  general  case ,  the  choice  of  tq  —  h/3/4  would  seem  to  be  a  reasonable,  or  at  least  an 
expedient,  choice  for  ro  for  more  complicated  unsymmetrical,  many-dimensional  systems. 

Before  preceeding  to  the  reactive  flux  formulation  based  on  eq.  (9) ,  some  consideration 
should  be  given  to  actual  evaluation  of  the  SCTST  rate  constant  for  a  general  condensed 
phase  activated  rate  process.  Clearly,  the  direct  numerical  evaluation  of  the  phase  space 
integral  in  eq.  (9)  is  prohibitively  time-consuming  for  a  many-dimensional  system.  For¬ 
tunately,  the  standard  techniques  of  computer  simulation  in  statistical  mechanics  can  be 
used.  For  example  one  can  formulate  the  evaluation  of  eq.  (9)  in  terms  of  a  Monte  Carlo 
average  by  using  a  reference  Hamiltonian,  HTef(r).  Equation  (9)  is  then  rewritten  as 

[o  rh.p/2—ro 

— —  J  drAif(r) 

where  the  average  ■  -)ref  denotes  averaging  in  the  semiclassical  flux- weighted  distribution 
function  for  the  reference  system  [cf.  eq.  (16)  below]  and  A H(t)  in  eq.  (12)  is  given  by 


ksCTST  —  k SCTST 


A H{t)  =  H\p(r ) , g(r), P(r) , Q(r )]  -  i7re/[pre/(r),gre/(r),Pre/(r),Qre/(r)]  .  (13) 


t 

\ 


i 


The  subscripted  variables  here  refer  to  the  reference  system  Hamiltonian  and  the  trajec¬ 
tories  on  that  inverted  potential  [cf.  eq.  (4)  written  for  the  reference  system],  while  the 
unsubscripted  variables  refer  to  the  actual  system.  Equation  (12)  can.  of  course,  be  written 
as  a  product  over  a  series  of  intermediate  reference  states  in  the  form  of  a  thermodynamic 
charging  algorithm.  As  an  alternative  to  the  reference  state/charging  procedure,  one  could 
instead  employ  umbrella  sampling  [16]  to  calculate  the  flux-weighted  distribution  function 
relative  to  the  reactant  state. 

As  demonstrated  by  Chapman.  Garrett,  and  Miller  [8],  the  SCTST  formula  in  eq.  (9) 
can  describe  the  tunneling  effects  in  gas  phase  atom-diatom  activated  rate  processes  and, 
as  discussed  in  the  next  section,  there  is  every  reason  to  expect  the  same  will  be  true  in  the 
condensed  phase.  Chapman,  et  al.  also  mentioned  that  one  could  run  classical  trajectories 
to  estimate  the  dynamical  recrossing  effects.  These  trajectories  would  be  initiated  from  the 
top  of  the  barrier  with  initial  conditions  weighted  by  the  semiclassical  distribution  function 
in  eq.  (7).  Although  CGM  did  not  carry  out  such  a  calculation,  we  believe  this  suggestion 
has  considerable  merit,  particularly  in  light  of  the  more  recent  development  of  the  classical 
reactive  flux  correlation  function  formalism  [3].  The  latter  formalism,  which  readily  carries 
over  to  the  semiclassical  case,  allows  one  to  efficiently  calculate  the  dynamical  corrections 
to  the  classical  TST  rate  constant. 

In  order  to  develop  a  semiclassical  reactive  flux  correlation  function  expression,  one 
must  first  define  a  dynamical  correction  factor  ksc  to  the  SCTST  rate  constant  in  eq.  (9), 
i.e., 

ksc  =  k  sc  k  SCTST  ,  (14) 

where  ksc  is  a  semiclassical  forward  rate  constant  which  includes  the  transition  state 
recrossing  effects  relative  to  the  SCTST  value  of  the  rate  constant.  These  dynamical 
effects  are  quantified  in  the  value  of  the  semiclassical  dynamical  correction  factor  ksc , 
given  as  a  function  of  time  by 

*sc(t)  =  {h[q(t)  -  q*])j+c  -  (h[q(t)  -  q*])j7c  ,  (15) 
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where  the  symbol  denotes  classical  trajectories  \p(t),q(t),  P(t),  Q(f)]  which  have 

positive  (“+”)  or  negative  (“— ”)  initial  reaction  coordinate  momenta.  The  reactive  flux 
trajectories  in  eq.  (15)  have  initial  conditions  generated  from  a  semiclassical  flux-weighted 
distribution  function  Wfc(p,  q,  P,  Q),  given  by 

!nl  f  2  /•fc/V 2— To  ) 

tVfc(p,g,P,Q)  =  (QRh)-xl(9- -g)^expi~- drff(r)  ^  .  (16) 

It  should  be  noted  that  the  classical  trace  of  the  semiclassical  flux-weighted  distribution 
function  yields  the  SCTST  rate  constant  in  eq.  (9),  i.e., 

kscTST  =  rrcdttf^g.P.Q)]  .  (IT) 

The  semiclassical  reactive  flux  correlation  function  in  eq.  (15)  is  calculated  (i.e.,  the 
trajectories  are  propagated)  until  it  reaches  a  “plateau”  value  [3].  After  that  point,  the 
value  of  ksc  is  essentially  constant  and  can  be  easily  determined.  It  should  also  be  noted 
that  the  dividing  surface  in  eq.  (16)  can  be  optimized  to  minimize  the  recrossing  effects 
and  hence  improve  the  accuracy  of  the  SCTST  rate  constant  in  eq.  (17).  However,  the 
semiclassical  flux  through  different  dividing  surfaces  may  not  be  invariant  with  respect  to 
the  choice  of  the  dividing  surface.  This  inconsistency  can  be  readily  eliminated,  albeit  at 
some  cost  in  computational  effort  (see  the  Appendix). 

Equations  (14)  and  (15)  represent  a  unified  computational  methodology  for  describing 
tunneling  effects  (via  k sctst)  and  dynamical  recrossing  effects  (via  ksc )  in  thermally 
activated  rate  processes.  There  is  no  doubt  that  much  of  the  CGM  formulation,  as  well 
the  present  reactive  flux  extension  of  that  theory,  has  a  somewhat  ad  hoc  character  (i.e.,  it 
has  not  been  derived  from  first  principles).  To  formulate  the  above  equations,  or  equations 
like  them,  using  a  completely  systematic  semiclassical  analysis  will  require  a  new  and  better 
understanding  of  the  quantum-classical  correspondence  principle.  This  rather  substantial 
challenge  will  be  reserved  for  future  research.  However,  as  will  be  shown  in  the  next  section, 
the  above  semiclassical  methodology  also  has  many  appealing  characteristics  within  the 
context  of  quantum  activated  dynamics,  not  the  least  of  which  is  its  clear  computational 
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power  for  studying  highly  nonlinear,  complex  systems.  In  fact,  it  is  the  latter  feature  of 
the  method  which  mitigates  to  some  degree  the  lack  of  rigor  in  its  formulation. 

3.  Applications 

3.1.  One- dimensional  Barrier  Crossing  Problems 

Much  about  the  characteristics  of  the  basic  SCTST  formula  in  eq.  (9)  can  be  un¬ 
derstood  by  analysis  of  a  simple  one-dimensional  barrier  crossing  problem.  For  example, 
CGM  have  already  demonstrated  that  eq.  (9),  for  the  choice  of  tq  =  ft/3/ 4,  gives  the  exact 
result  when  the  barrier  potential  14 (q)  is  approximated  to  be  an  upside-down  parabola 
[i.e.,  14(g)  ~  Vo  —  (l/2)mu >l(q  —  g*)2].  This  result,  however,  is  essentially  a  “higher” 
temperature  result,  being  valid  only  for  temperatures  such  that  h,(3u>b  <  2ir. 

For  temperatures  below  the  one-dimensional  crossover  temperature  (i.e.,  for  h/3ub  > 
27t),  the  SCTST  formula  in  eq.  (9)  can  be  tested  numerically.  Indeed,  CGM  have  presented 
results  for  the  Eckart  barrier  and  for  the  coilinear  H  ~f  H2  reaction  which  were  in  good 
agreement  with  the  known  exact  results.  But  why  the  good  agreement?  In  this  regard, 
it  proves  to  be  quite  instructive  to  carry  out  a  more  in-depth  analysis  of  the  dominant 
contributions  to  the  phase  space  integral  in  eq.  (9)  for  a  non-quadratic  barrier  potential 
such  as  a  symmetric  Eckart  barrier  [14(g)  =  Vocosh~2(atq/2)}.  For  this  (or  any)  one¬ 
dimensional  system,  the  integral  over  the  reaction  coordinate  q  in  eq.  (9)  is,  of  course, 
fixed  at  q*  by  the  delta-function.  The  integral  over  the  momentum  must  be  performed, 
however.  These  values  of  momenta  provide  the  initial  conditions  for  the  trajectories  which 
are  propagated  on  the  upside-down  barrier  potential,  —14(g),  in  the  calculation  of  the 
exponential  weighting  function  in  eq.  (9).  In  fig.  1,  the  flux- weighted  distribution  function 
is  plotted  for  a  symmetric  Eckart  barrier  of  height  2000  cm-1  and  lengthscale  a  =  5.73 
A"1  at  125  K.  It  is  seen  that  the  distribution  function  is  strongly  peaked  at  a  certain 
value  of  reaction  coordinate  momentum  given  by  pmai.  In  fig.  2,  the  imaginary  time 
trajectory  [cf.  eq.  (4)]  on  the  upside-down  Eckart  barrier  potential  is  plotted  for  the  initial 
conditions  (pmax  ><?*)>  where  q*  equals  zero  in  this  case.  From  this  initial  condition,  the 
trajectory  is  propagated  backward  in  time  to  the  time  —  to  =  —ft/3/4  and  forward  in  time 
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to  ro  =  TiQ/ 4.  It  is  seen  from  fig.  2  that  the  trajectory  exactly  reaches  the  classical  turning 
point  (dashed  lines)  of  the  inverted  barrier  potential.  Thus,  this  trajectory  is  precisely  the 
WKB  trajectory  for  tunneling  through  the  barrier.  Stated  differently,  the  trajectory  shown 
in  fig.  2  is  exactly  one-half  of  the  well-known  semiclassical  periodic  orbit,  or  “kink-pair”, 
trajectory  from  instanton  theory  [7,9].  (It  should  be  noted  that  the  SCTST  expression  in 
eq.  f°)  requires  only  one- half  of  an  instanton  periodic  orbit  because  of  the  factor  of  two  in 
the  exponent.) 

Actually,  it  is  no  accident  that  the  instanton  dynamics  dominates  the  value  of  the 
SCTST  formula  in  eq.  (9)  at  low  temperatures.  In  general,  a  steepest  descents  analysis 
of  the  exponential  in  eq.  (9)  below  the  crossover  temperature  will  locate  two  dominant 
contributions  to  that  term:  the  aforementioned  instanton  trajectory  and  the  “trivial” 
stationary  trajectory  solution  for  the  initial  conditions  p  =  0,  q  =  q* .  However,  the  latter 
contribution  to  the  exponential,  which  is  quite  large,  cannot  contribute  to  the  flux- weighted 
distribution  function  in  eq.  (16)  which  underlies  the  SCTST  expression  [eq.  (9)].  This 
feature  of  the  theory  arises  because  the  term  [pf  in  the  flux-weighted  distribution  insures  a 
weight  of  zero  for  the  p  =  0,  q  =  q*  contribution  to  the  semiclassical  phase  space  integral. 
As  a  result  of  this  feature,  the  basic  SCTST  formula  in  eq.  (9)  will  always  have  a  value 
very  close  to  the  instanton  result  for  temperatures  below  crossover.  Thus,  in  view  of  this 
low  temperature  analysis  and  the  higher  temperature  parabolic  barrier  analysis  [8],  it  can 
be  concluded  that  the  SCTST  formula  can  be  faithfully  used  at  all  temperatures ,  providing 
a  stable,  continuous  theory  for  the  rate  constant  above,  below,  and  through  the  crossover 
region.  (This  very  appealing  property  is  also  a  characteristic  the  theory  has  in  common 
with  the  path  integral  QTST  [10]).  In  order  to  highlight  the  accuracy  of  the  SCTST 
formula  in  eq.  (9)  as  a  function  of  temperature,  Table  I  lists  results  computed  with  that 
formula  for  the  symmetric  Eckart  barrier  quantum  corrections  versus  the  known  exact 
results  [17]. 

3.2.  A  Two-Dimensional  Test  of  the  Reactive  Flux  Method 

The  previous  subsection  contains  a  discussion  of  the  ability  of  the  basic  SCTST  for¬ 
mula  in  eq.  (9)  to  accurately  describe  tunneling  effects.  In  this  subsection,  a  test  of  the 
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semiclassical  reactive  flux  formulation  in  eqs.  (14)— (16)  is  presented.  In  order  to  do  this, 
the  reactive  flux  was  studied  for  a  non-quadratic  barrier  potential  which  is  bilinearly  cou¬ 
pled  to  a  harmonic  oscillator.  Exact  quantum  reactive  flux  results  exist  for  this  system 
[18],  and  it  has  also  been  studied  by  other  approximate  methods.  The  Hamiltonian  used 
in  these  studies  is  given  by 

H  =  +  +  W  +  -  c«)2  ,  (18) 

where  p  is  a  reduced  mass,  c  is  the  coupling  between  modes,  and  u  the  bath  mode  frequency. 
The  reaction  coordinate  potential  was  chosen  to  be  a  symmetric  Eckart  potential 

Vb(q)  =  V0sech2(aq/2)  .  (19) 

where  a  is  a  length  parameter,  VQ  is  the  barrier  height,  and  the  frequency  of  motion  on 
the  top  of  the  Eckart  barrier  is  defined  by  u2  =  —a2  Vo/ (2/z). 

The  coupling  constant  c  in  eq.  (18)  was  determined  by  the  first  term  in  a  cosine 
Fourier  representation  of  a  Gaussian  friction  kernel  such  that  [18] 

^  2  rTbath 

c 2  = — r -  /  dtcos(ut)r](t ) 

^  Tbath  Jo 

_  2/M O'2  exp[-{ua)2/2]  ^ 

Tbath  (wcr)2 

where  c  determines  the  timescale  of  the  friction  kernel.  The  variable  /  here  is  a  reduced 
friction  parameter,  given  by  T]o/p\u>b\,  and  rjo  is  the  zero  time  value  of  the  friction. 

In  the  calculations,  the  barrier  height  was  taken  to  be  3425.3  cm-1,  the  lengthscale  a 
to  be  3.97  A-1,  the  reduced  mass  p  equal  to  0.672  amu,  and  the  barrier  frequency  u>b  to  be 
1164.0  cm-1.  For  the  bath,  the  decay  parameter  a  was  chosen  to  be  86.6  fs  and  the  bath 
fequency  u  to  be  48  cm-1.  The  timescale  parameter  Tbath  in  eq.  (20)  was  set  equal  to  2a. 
(This  model  is  referred  to  as  the  “slow  friction”  in  ref.  [18].)  In  the  calculations,  2000  sets 
of  initial  conditions  were  generated  using  the  SC  flux  weighted  distribution  function  in  eq. 
(16).  From  this  set  of  initial  conditions,  2000  trajectories  were  then  run,  each  for  1  ps,  to 
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calculate  use  from  eq.  (15).  The  calculations  were  performed  at  temperatures  of  200K 
and  300K. 

As  can  be  seen  from  Figs.  3  and  4,  the  combined  semiclassical  theory  in  eqs.  (14)— (16) 
(i.e.,  SCTST  with  reactive  flux  corrections)  is  in  very  good  agreement  with  the  exact 
quantum  results.  On  the  other  hand,  the  SCTST  rate  constant  alone  is  not  in  good 
agreement  with  the  exact  results.  One  can  conclude,  therefore,  that  the  inclusion  of  the 
semiclassical  dynamical  correction  factor  in  eq.  (15)  is  essential  for  the  accuracy  of  the 
method  at  large  couplings.  At  intermediate  couplings,  the  semiclassical  reactive  flux  result 
drops  drops  too  far  below  the  exact  results,  particularly  at  200  K.  This  regime  is  due  to  an 
interesting,  but  pathological,  feature  of  this  two-mode  problem  which  leads  to  metastable 
trapped  classical  trajectories  in  the  transition  state  region.  Such  trapping  regions  in  phase 
space  are  far  more  pronounced  in  two-dimensional  potentials  for  which  locally  integrable 
classical  motion  can  rigorously  segment  phase  space  into  disjoint  regions  [19].  In  many¬ 
dimensional  condensed  phase  systems,  this  behavior  will  be  far  less  likely  to  occur. 

3.3.  Parabolic  Barrier  Coupled  to  a  Multidimensional  Bath  Model 

To  this  point  in  the  present  manuscript,  the  SCTST  rate  constant  and  SC  reactive 
flux  method  have  been  applied  to  low-dimensional  test  cases.  These  models  were  chosen,  of 
course,  because  exact  results  are  essentially  non-existent  for  xnany-body  quantum  activated 
rate  problems.  One  exception  to  this  rule,  however,  is  the  many-dimensional  parabolic 
barrier  potential.  This  potential  is  given  by 

V(q,  Q)  =  V0  -  -rrwlq2  +  (  Qi  +  — —q)  ,  (21) 

2  2  V  miUi  / 

where  u n,  is  the  magnitude  of  the  unstable  barrier  frequency,  the  constants  c,  are  the 
coupling  constants  between  the  reaction  coordinate  and  the  bath  modes,  and  the  bath 
mode  frequencies  and  masses  are  given  by  u>x  and  mt,  respectively.  The  bath  modes  in  the 
above  Hamiltonian  are  chosen  to  model  the  linear  force  fluctuations  from  a  bath  on  the 
reaction  coordinate  [20]. 

For  relatively  few  degrees  of  freedon,  the  multidimensional  parabolic  barrier  problem 
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is  trivial  because  a  normal  coordinate  rotation  can  be  readily  performed  to  separate  the 
true  unstable  reactive  normal  mode  p  from  the  stable  orthogonal  vibrations  {y} }.  However, 
for  an  infinite  number  of  bath  oscillators  as  in  condensed  phase  problems,  the  normal  mode 
transformation  can  only  be  defined  formally.  Even  so,  if  such  a  transformation  is  performed 
to  find  the  many-dimensional  set  of  normal  modes,  then  the  SCTST  in  eq.  (9)  can  still  be 
applied  to  the  Hamiltonian  arising  from  eq.  (21),  but  with  the  dividing  surface  chosen  to 
be  along  the  true  multidimensional  reactive  mode  p  [21,22]  rather  than  along  the  “simple” 
reaction  coordinate  q  as  is  written  in  eq.  (9)  In  this  case,  the  result  for  the  rate  constant 
can  be  obtained  analytically  and  is  given  by  the  exact  result  [22]: 


ksCTST 


kBT 


Wj/2 


N0 


n 


exp(—/3Vo) 


(22) 


hQ*  sin(hfi\ll 2)  fJi  2s inh(hpXj/2 ) 

where  Aq  is  the  magnitude  of  the  frequency  of  the  multidimensional  reactive  normal  mode 
and  the  frequencies  {Aj}  are  the  frequencies  of  the  stable  modes  orthogonal  to  the  reaction 
coordinate.  The  continuum  (Nosc  — ►  oo)  limit  of  the  multidimensional  parabolic  result  is 
given  by  [22,23] 


*SCTST  =  ££  'I  II 


exp(-/3Vo ) 


hQn  ^  -  <4  +  SM(SWm  ' ' ' " ' ' "  ’  (23) 

where  Cln  equals  27rn/fi/3,  t)(-  •  •)  is  the  Laplace  transform  of  the  friction  kernel  along  the 
coordinate  q,  and  Aq  is  given  from  the  solution  of  the  well-known  Grote-Hynes  relationship 
[24] 


4  =  -7— 1 —t -  •  (24) 

The  result  in  eq.  (23)  can  be  shown  [23]  to  provide  an  accurate  description  of  quantum 
activated  dynamics  for  reactive  systems  coupled  to  linearly  responding  condensed  phase 
environments  which  are  above  the  multidimensional  crossover  temperature  (i.e.,  for  Ti/3Xq  < 
27r). 

The  result  in  eq.  (23)  is,  of  course,  the  exact  result  [22,23].  It  is  therefore  concluded 
that  a  rotation  of  the  dividing  surface  in  the  basic  SCTST  rate  constant  expression  can 
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improve  the  accuracy  of  the  theory-at  least  for  the  parabolic  sy.-tem.  The  question  arises, 
however,  as  to  whether  a  calculation  using  the  semiclassical  reactive  flux  method  with  a 
flux-weighted  distribution  function  for  an  unrotated  dividing  surface  [cf.  eq.  (16)]  will 
achieve  the  same  result.  As  it  turns  out,  such  a  calculation  will  always  give  the  exact  rate 
constant.  To  prove  this,  one  first  transforms  the  multidimensional  quadratic  Hamiltonian 
for  the  potential  in  eq.  (21)  into  the  separable  normal  mode  form  [21,22] 


where,  for  convenience,  Vo  in  eq.  (21)  is  taken  to  be  zero.  Since  the  Hamiltonian  is 
separable  in  this  coordinate  system,  one  can  readily  uncover  the  initial  conditions  for  the 
trajectories  which  will  always  proceed  from  reactants  to  products.  These  conditions  are 
represented  by  the  characteristic  function  x(Pp,  p)  —  h(pP)h(p2p- Xq2  p2) ,  where  h(-  ■  •)  is  the 
unit  step  function.  This  function  simply  states  that  any  trajectory  with  an  initial  positive 
momentum  along  the  true  multidimensional  reaction  coordinate  p  and  an  energy  above 
the  barrier  along  that  coordinate  will  be  reactive  (i.e.,  it  will  contribute  a  weight  of  unity). 
The  calculation  of  the  semiclassical  reaction  rate  from  eqs.  (14)  —  (17)  is  then  equivalent 
to  averaging  the  characteristic  function  with  the  semiclassical  flux-weighted  distribution 
function  in  eq.  (16).  A  somewhat  lengthy,  but  straightforward,  evaluation  of  Gaussian 
integrals  then  leads  to  the  result  that  the  semiclassical  reactive  flux  calculated  with  the 
semiclassical  procedure  specified  by  eqs.  (14)  —  (17)  will  be  exact  for  the  many-dimensional 
parabolic  model. 

Results  are  shown  in  Fig.  5  for  the  numerically  determined  (exact)  transmission 
coefficient  for  the  parabolic  barrier  model  with  a  barrier  frequency  of  1000  cm-1  and  a  100 
discrete  oscillator  bath  [25].  The  bath  was  chosen  to  model  an  exponential  friction  kernel 
with  a  decay  constant  of  4u;^1.  The  strength  of  the  friction  is  characterized  here  by  the 
dimensionless  quantity  r)(0)/mujl,  where  tj(0)  is  the  zero  time  value  of  the  friction  kernel. 
The  reactant  partition  function  was  assumed  to  be  given  by  a  harmonic  oscillator  potential 
along  the  q  direction  with  a  frequency  of  V2u>b  which  is  also  coupled  to  the  oscillator  bath. 
The  numerical  stability  and  relative  simplicity  of  this  multidimensional  calculation  bodes 
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well  for  applications  of  the  semiclassical  method  to  other,  more  complicated,  many-body 
systems— at  least  above  the  crossover  temperature. 


3.4.  Escape  from  a  Metastable  Well  in  a  Dissipative  Medium 

Rips  and  Poliak  [26]  have  presented  a  theory  for  the  quantum  mechanical  escape 
rate  of  a  particle  out  a  metastable  well  in  the  presence  of  a  dissipative  medium.  Their 
theory  is  valid  in  the  case  of  a  linearly  responding  bath  and  for  temperatures  above  the 
multidimensional  crossover  temperature.  The  Rips-Pollak  expression  for  the  quantum 
mechanical  escape  rate  is  given  by  [26] 

t  =  Texpt-m  ,  (26) 

where  Vo  is  the  barrier  height,  ujq  is  the  well  frequency,  u>b  is  the  barrier  frequency,  and  Aq 
is  the  Grote-Hynes  reactive  frequency  [cf.  eq.  (24)].  The  factor  E  in  eq.  (26)  accounts  for 
the  particle  tunneling  in  the  presence  of  the  dissipation  and,  following  the  notation  in  eq. 
(23) ,  is  given  by 


n 


+  u>0  + 


LJi  -  U)l  -f  nnf)(On)/m 


(27) 


The  important  new  feature  in  the  Rips-Pollak  expression  is  the  quantum  mechanical  de¬ 
population  factor  T,  given  by 


T  =  exp 


0hxlsin((3hxl/2) 

2t r 
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ln[  1  —  P(t  —  i/ 2)] 

00  cosh(rf3h\o )  -  cos(j3fi Aq/2) J 


dr 


(28) 


In  this  equation,  P(r  —  i/2 )  is  the  Fourier  transform  of  the  quantum  mechanical  transition 
probability,  given  in  terms  of  the  set  of  stable  normal  mode  frequencies  {Aj}  by 


P(t  —  i/2 )  =  exp 


No 


-E- 


Vi[cosh(l3h\i/2)  -  cos(T/3h\i )] 


»= 1 


smh(0JhXi/2 ) 


(29) 


where  the  quantities  ut  are  given  by 
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Vt 


dt\etXitl  F(t\) 


(30) 


The  term  F(t )  here  is  the  time-dependent  zero-order  force  acting  on  the  unstable  normal 
mode  p  such  that 


F(t)  =  -  uoom~l/2 


dVNL(q) 


dq 


q=m~1/3uoop 


(31) 


where  V)vL(g)  is  the  nonlinear  part  of  the  barrier  potential,  defined  as  Vnl  —  Vb(q)  + 
(l/2)nujlq‘2 .  The  projection  of  the  barrier  coordinate  q  onto  the  unstable  multidimensional 
normal  mode  p  and  the  stable  normal  modes  {t/,}  is  determined  by  the  transformation  [21] 


_  N0.c 

y/rnq  =  uqqP  +  ^  •  (32) 

i=l 

The  constants  gt  in  eq.  (30)  are  defined  by  the  relation  =  u,o/«oo-  More  details  about 
the  Rips-Pollak  theory  can  be  found  in  ref.  [26]. 

The  present  semiclassical  reactive  flux  formalism  was  applied  to  a  cubic  barrier  po¬ 
tential  coupled  to  a  linearly  responding  bath  with  an  exponential  friction  kernel.  The 
potential  along  q  was  given  by 

Vb(q)  =  Vo-  ^ mwlq 2  -  | q 3  ,  (33) 

where  a2  =  (ma;^)3/6Vo.  In  the  simulations,  the  barrier  height  Vo  was  given  a  value  of 
3425.3  cm-1  and  the  relevant  particle  frequencies  were  chosen  to  be  Uf,  =  uo  =  1000  cm-1. 
The  linearly  responding  bath  was  modeled  by  100  harmonic  oscillators  as  in  our  earlier 
work  on  classical  activated  dynamics  [25].  The  exponential  decay  timescale  of  the  bath 
friction  kernel  was  chosen  to  be  4  u;^1.  As  in  Sec.  3.3,  the  strength  of  the  friction  was 
characterized  by  the  dimensionless  quantity  77(0) /mw2,  where  77(0)  is  the  zero  time  value 
of  the  friction  kernel. 

The  results  from  the  semiclassical  reactive  flux  calculations  for  the  cubic  well  model  are 
depicted  in  Fig.  6  and  compared  to  the  analytic  results  from  the  Rips-Pollak  theory.  Also 
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shown  for  comparison  are  the  analytic  classical  result  from  the  theory  of  Poliak,  Grabert, 
and  Hanggi  (PGH)  [27].  As  in  the  classical  theory,  there  are  two  mechanistic  pictures 
required  to  analyze  the  escape  rate  of  the  particle  from  the  well  as  a  function  of  friction 
strength  [2].  At  low  friction,  the  dynamics  is  determined  by  the  rate  of  energy  exchange 
between  the  particle  in  the  well  and  the  bath.  As  the  coupling  increases  in  this  regime,  the 
rate  increases.  At  high  frictions,  thermal  equilibrium  is  maintained  in  the  well  and  the  rate 
is  determined  by  passage  over  (or  through)  the  barrier.  Therefore,  at  high  frictions,  the 
rate  decreases  as  the  coupling  increases.  This  phenomenology  is,  of  course,  the  standard 
Kramers  picture  of  activated  dynamics  [28],  and  the  quantum  mechanical  modifications  of 
that  behavior  are  evident  from  Fig.  6.  The  agreement  between  the  SC  reactive  flux  and 
Rips-Pollak  theory  is  quite  encouraging  and  again  suggests  that  the  semiclassical  theory 
will  be  useful  for  complicated  condensed  phase  problems. 

3. 5.  Double  Well  in  a  Linearly  Responding  Dissipative  Bath 

As  a  final  test  of  the  SCTST  and  the  semiclassical  reactive  flux  formalism,  the  rate 
constant  was  calculated  in  the  region  of  the  crossover  temperature  for  a  double  well  barrier 
potential  coupled  to  a  linearly  responding  bath.  The  crossover  temperature  in  this  case 
can  be  estimated  by  the  usual  relationship 

Tcr  =  h\l/(2irkB)  , 

where  Aq  is  the  Grote- Hynes  reactive  frequency  in  eq.  (24). 
double  well  was  employed  of  the  form 

Vb(q)  =  Vo-  q4  .  (35) 

This  quartic  double  well  has  the  property  that  the  well  frequency,  u;o,  is  related  to  the 
barrier  frequency,  Ub,  such  that  uo  =  y/2vb- 

The  specific  values  of  the  paramters  were  again  u>b  =  1000  cm-1  and  Vo  =  3425.35 
cm-1.  The  bath  was  again  modeled  by  100  discrete  oscillators  [25],  having  a  friction  kernel 
with  timescale  4u^J  and  strength  r/o//^  =  1.0.  For  this  value  of  the  friction,  the  unsta- 


(34) 


In  these  calculations,  a  quartic 
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ble  normal  mode  frequency  Ag  is  ~  530  cm'1;  the  corresponding  crossover  temperature 
being  ~  125  K.  Semiclassical  calculations  were  then  performed  at  100K.  125K,  and  150K. 
The  results  for  the  semiclassical  rate  constant  [eq.  (14)]  and  the  semiclassical  transmission 
coefficient  [eq.  (15)]  are  tabulated  in  Table  II.  For  the  sake  of  comparison,  the  semiclassical 
transmission  coefficient  at  300  K  was  found  to  have  a  value  of  0.68.  There  are  no  exact 
quantum  results  here  with  which  to  compare  the  semiclassical  calculations,  but  the  semi¬ 
classical  transmission  coefficient  exhibits  an  interesting  drop  in  value  as  the  temperature 
is  lowered  through  the  crossover  temperature.  This  result  is  in  general  agreement  with 
the  behavior  predicted  by  an  analytic  theory  [13]  for  the  transmission  coefficient  in  path 
integral  QTST  below  the  crossover  temperature. 


4.  Concluding  Remarks 

In  the  present  paper,  a  semiclassical  transition  state  theory  proposed  by  Chapman, 
Garrett,  and  Miller  [8]  has  been  “rediscovered”  and  analyzed  within  the  context  of  con¬ 
densed  phase  reaction  rate  theory.  For  the  cases  studied,  it  was  found  that  the  basic  SCTST 
formula  [eq.  (9)]  is  capable  of  accurately  treating  activated  rate  processes  over  a  wide  range 
of  temperatures,  including  in  the  strongly  tunneling  regime.  The  relationship  of  the  theory 
to  the  more  widely  known  instanton  theory  [7,9]  has  also  been  elaborated.  A  significant 
extension  of  the  theory  to  include  dynamical  dividing  surface  recrossing  effects  has  been 
presented  in  the  form  of  a  semiclassical  reactive  flux  correlation  function  formalism  [eqs. 
(14)— (16)].  The  accuracy  of  this  formalism  was  tested  against  several  model  and  analytical 
results  and  found  to  be  quite  good.  The  great  strength  of  the  method  is  the  ease  with 
which  it  can  be  numerically  applied  to  nonlinear,  many-body  systems.  Moreover,  many  of 
the  molecular  dynamics  and  Monte  Carlo  sampling  techniques  from  statistical  mechanics 
can  be  combined  with  the  formalism.  The  method  is  likely  to  be  broadly  applicable  to  a 
variety  of  important  problems,  including  proton  transfer,  adiabatic  electron  transfer,  and 
the  diffusion  hydrogen  on  material  surfaces.  Some  of  these  applications  will  be  reported 
in  future  publications. 
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Appendix 

As  stated  in  the  main  text,  one  can  imagine  an  optimization  of  the  dividing  surface  in 
eq.  (16)  to  minimize  recrossing  effects  in  the  calculation  of  the  semiclassical  transmission 
coefficient  [eq.  (15)].  In  such  a  procedure,  however,  the  semiclassical  flux  may  not  be 
invariant  with  respect  to  the  choice  of  the  dividing  surface— as  it  must  be  in  an  exact 
calculation.  Although  this  negative  aspect  of  the  formulation  will  probably  not  introduce 
major  errors,  one  can  in  fact  modify  the  method  to  make  the  flux  independent  of  the 
dividing  surface.  In  order  to  do  this,  a  semiclassical  effective  Hamiltonian  can  be  defined 
as 


ffts£  =  -fcBrfn[W-sc(p,q)]  .  (-41) 

where  the  semiclassical  density  Wscte,  q)  is  given  by  eq.  (7).  This  effective  Hamiltonian 
is  then  used  to  generate  the  recrossing  trajectories  in  the  calculation  of  the  semiclassical 
transmission  coefficient  [eq.  (15)].  The  equations  of  motion  are  given  in  this  case  by 
Hamilton’s  equations,  but  with  the  classical  Hamiltonian  replaced  by  the  semiclassical 
effective  Hamiltonian  H^ff-  Clearly,  solving  these  equations  will  be  more  computationally 
demanding  than  the  usual  classical  equations. 
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TABLE  I.  Comparison  of  Tunneling  Corrections  for  the  Symmetric  Eekart  Barrier. 


Temperature 

SCTST  Results* 

Exact  Results6 

125.6 

1686.0 

1970.0 

150.7 

144.4 

162.0 

188.3 

20.3 

22.0 

251.1 

4.9 

4.2 

301.3 

2.9 

3.1 

376.7 

2.0 

2.1 

502.2 

1.5 

1.5 

753.3 

1.2 

1.2 

“Equation  (9) 
kFrom  ref.  [17]. 
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TABLE  II.  Rate  Results  for  the  Double  Well  in  a  Dissipative  Bath. 


Temperature 

«sca 

Rate  (s  l)b 

150 

0.91 

1.5  x  10"2 

125 

0.87 

7.9  x  10“4 

100 

0.40 

2.0  x  10“8 

“Equation  (15). 

The  numerical  error 

is  5%  on  the  two  higher 

temperature  results  and  20%  on  the  lowest. 
bEquation  (14) 
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Figure  Captions 


Fig.  1:  A  plot  of  the  semiclassical  flux-weighted  distribution  function  [eq.  (16)]  as  a  function 
of  initial  momentum  for  a  symmetric  Eckart  barrier  at  T  =  125 K. 

Fig.  2:  The  imaginary  time  trajectory  on  the  upside  down  Eckart  potential  corresponding 
to  the  initial  condition  (pmai ,  0)  which  give  the  maximum  value  of  the  flux- weighted 
distribution  function  in  Fig.  1.  The  trajectory  is  propagated  forward  in  time  to  ft/3/4 
and  backward  in  time  to  —ft/3/4.  The  dashed  lines  are  the  classical  turning  points  on 
the  inverted  Eckart  barrier  potential  for  this  initial  condition. 

Fig.  3:  Plots  of  the  reactive  flux  for  the  two-dimensional  model  in  Sec.  3.2.  The  results 
are  from  the  SCTST  expression  [eq.  (9)],  the  semiclassical  expression  which  includes 
the  semiclassical  transmission  coefficient  [eq.  (14)],  and  exact  quantum  wavepacket 
calculations  [18]  at  200K. 

Fig.  4:  Same  as  in  Fig.  3,  but  for  a  temperature  of  300K. 

Fig.  5:  A  plot  of  the  numerically  determined  semiclassical  rate  constant  [eq.  (14)]  for  a 
parabolic  barrier  potential  linearly  coupled  to  100  bath  oscillators  as  a  function  of  the 
friction.  The  solid  line  is  the  exact  analytic  result.  The  reactant  partition  function  is 
assumed  here  to  be  for  a  harmonic  oscillator  potential  coupled  to  the  same  bath  as  at 
the  barrier. 

Fig.  6:  A  comparison  of  the  rate  constant  for  a  metastable  cubic  well  in  a  dissipative  100 
oscillator  bath  as  calculated  from  the  Rips-Pollak,  PGH,  and  semiclassical  [eq.  (14)] 
theories  at  300K  (see  Sec.  3.4). 


-25- 


OJx]  (x‘b‘d)jj. 


[*n*«]  (x)b 


IO 

C 


CO 


CM 


[i-s]  xnu  8AIJ0B8H 


300K 


r-  co  10  co 


[i-S]  xny  9Al^0B8^ 


i a 


CM 


lA 


[i-s]  9^a 


CM 


